
function [icell,jcell,xcorn,ycorn] = xy2ij_curv (x,y,x_grid,y_grid,isearch) ;
%
%-----------------------------------------------------------------------------
%
% 	Finds model grid cell indices (i,j) belonging to a set of (x,y) points
%	for a (very) curvilinear grid.
%
%-----------------------------------------------------------------------------
%
%   Usage :
%	    [icell,jcell,xcorn,ycorn] = xy2ij_curv (x,y,x_grid,y_grid,isearch)
%
%   Input :
%	    x = column vector containing x coords. of points of interest.
%	    y = column vector containing y coords. of points of interest.
%	    x_grid = matrix of dimension (ncelly+1,ncellx+1) containing
%                   the x coords. of the grid points.
%	    y_grid = ditto for y.
%	    isearch = radius in grid cells to do cell search from nearest
%		      grid point.
%
%   Output :
%	    icell = column vector containing the cell i values (starting at 1).
%	    jcell = ditto for j.
%           (NaN,NaN) is returned if the point lies outside the grid.
%	    xcorn = matrix of dimension (4,npt) of x coords. of cell corners.
%	    ycorn = matrix of dimension (4,npt) of y coords. of cell corners.
%
% Note that this function is part of PLUM (Plot Lots Using Matlab) written  
% by John Andrewartha 
% Author: John Andrewartha , CSIRO
% email: John.Andrewartha@csiro.au
% Feb 2013 
%-----------------------------------------------------------------------------
%

	[m,n] = size(x_grid) ;
	[npt] = length(x) ;
%
	icell = NaN * ones(npt,1) ; jcell = NaN * ones(npt,1) ;
	xcorn = NaN * ones(4,npt) ; ycorn = NaN * ones(4,npt) ;
%
% LOOP OVER THE POINTS
%
	for ip = 1:npt,
	    xg = x_grid-x(ip) ; yg = y_grid-y(ip) ;
	    dist2 = xg.^2 + yg.^2 ;
	    [ymin,jmin] = min(dist2) ;
	    [dmin,i] = min(ymin) ;
	    j = jmin(i) ;
%
	    imin = max([1 i-isearch]) ;
	    jmin = max([1 j-isearch]) ;
	    imax = min([n-1 i+isearch]) ;
	    jmax = min([m-1 j+isearch]) ;
%
	    for i = imin:imax,
		for j = jmin:jmax,
	            xpoly = [xg(j,i);xg(j,i+1);xg(j+1,i+1);xg(j+1,i);xg(j,i)] ;
	            ypoly = [yg(j,i);yg(j,i+1);yg(j+1,i+1);yg(j+1,i);yg(j,i)] ;
		    isin = 0 ;
		    if (~isnan(sum(xpoly)+sum(ypoly))),
%...................... Do a point-in-poly.
	   	        a = ypoly(1:4) ;
	    	        b = ypoly(2:5) ;
	                ind1 = find(a>=0 & b<0) ;
	                ind2 = find(a<0 & b>=0) ;
	                cx = [xpoly(ind1) ; xpoly(ind2)] ;
	                cy = [ypoly(ind1) ; ypoly(ind2)] ;
	                dx = [xpoly(ind1+1) ; xpoly(ind2+1)] ;
	                dy = [ypoly(ind1+1) ; ypoly(ind2+1)] ;
	                if (length(cx)>0),
	        	    e = cx - cy.*(dx-cx) ./ (dy-cy) ;
	        	    ind3 = find(e>0) ;
	        	    if (length(ind3)>0) ,
	            	        ncross = length(ind3) ;
	                        isin = rem(ncross,2) ;
	                    end
	                end
	                if (isin),
	                    icell(ip) = i ; jcell(ip) = j ;
	                    xcorn(:,ip) = xpoly(1:4)+x(ip) ; ycorn(:,ip) = ypoly(1:4)+y(ip) ;
		            break ;
		        end
		    end
	            if (isin), break ; end
	        end
	        if (isin), break ; end
	    end
%
	    if(rem(ip,500)==0),
		disp(['xy2ij : done ' int2str(ip) ' points']) ;
	    end
	end
